1 Setup

1.1 Libs

library(car) 
## Loading required package: carData
library("Rmisc")
## Loading required package: lattice
## Warning: package 'lattice' was built under R version 4.3.1
## Loading required package: plyr
## Warning: package 'plyr' was built under R version 4.3.1
library("ggpubr")
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.3.1
## 
## Attaching package: 'ggpubr'
## The following object is masked from 'package:plyr':
## 
##     mutate
library("cowplot")
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggpubr':
## 
##     get_legend
library("lme4")
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 4.3.1
library("phyloseq")
library("glmmTMB")
## Warning: package 'glmmTMB' was built under R version 4.3.1
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.6.0
## Current Matrix version is 1.6.1.1
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
library("emmeans")
## Warning: package 'emmeans' was built under R version 4.3.1
library(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
#install.packages("multcompView")
library(multcompView)
library("bbmle")
## Loading required package: stats4
library("DHARMa")
## This is DHARMa 0.4.6. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library("nlme")
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:lme4':
## 
##     lmList
library("bestNormalize")
## 
## Attaching package: 'bestNormalize'
## The following object is masked from 'package:MASS':
## 
##     boxcox

1.2 Data

setwd("~/nicolagk@hawaii.edu - Google Drive/My Drive/Mosquito_business/Mosquito_microbes_git/04.microbiome_analysis/04a.diversity")

ps.clean <- readRDS("../../02.process_asvs/ps.clean.trim.rds")
#ps.clean <- readRDS("../../02.process_asvs/ps.clean.trim.rare9200.rds")

samdf <- data.frame(ps.clean@sam_data)

##made below
#df.div <- read.csv("div.trim.csv",row.names=1)

2 Alpha diversity calculations (first run through)

##option to run rarefied stuff:
#ps.clean <- ps.rare

df.all <- data.frame(estimate_richness(ps.clean, split=TRUE, measures=c("Shannon","InvSimpson","Observed")))

df.all$sample_name <- rownames(df.all)

df.all.div <- merge(df.all,samdf,by="sample_name") #add sample data

#shannon diversity divided by species richness
df.all.div$even <- df.all.div$Shannon/(log(df.all.div$Observed))

df.div <- df.all.div

#adding most recent sample counts
sums <- data.frame(sample_sums(ps.clean))
colnames(sums) <- c("lib_size_trim")
sums$sample_name <- row.names(sums)

df.div <- merge(df.div,sums,by="sample_name")

#write.csv(df.div,"div.trim.csv")

3 Analyzing alpha metrics

3.1 All data

Very strange/interesting pattern where the water types that shouldn’t have a lot of diversity have a lot of OTUs, but not a lot of reads…

ggplot(df.div,aes(x=type,y=Observed,color=infusion))+
  geom_boxplot()+
  geom_jitter()

ggplot(df.div,aes(x=type,y=totsum,color=infusion))+
  geom_boxplot()+
  geom_jitter(position=position_jitterdodge())

3.2 Plotz

#just experimental types
df.div.exp <- subset(df.div,type=="Microbial Water"|type=="A.albopictus")

df.div.exp$inf.temp <- paste0(df.div.exp$infusion,df.div.exp$temperature)

df.div.exp$type <- sub("A.albopictus","Mosquitoes",df.div.exp$type)
df.div.exp$type <- sub("Microbial Water","Meso. water",df.div.exp$type)

df.div.exp$type <- factor(df.div.exp$type,levels=c("Mosquitoes","Meso. water"))

3.2.1 Richness

df.div.exp.se <- summarySE(df.div.exp,measurevar="Observed",groupvars=c("inf.temp","infusion","temperature","type"))

##dots & error bars
gg.rich <- ggplot(df.div.exp,aes(x=infusion,y=Observed,fill=inf.temp,shape=inf.temp))+
  geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
  geom_errorbar(data=df.div.exp.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
  geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
  #geom_boxplot(outlier.shape=NA,alpha=0.5)+
  facet_wrap(~type)+
  theme_cowplot()+
  ylab("ASV richness")+
  xlab("Infusion")+
  scale_x_discrete(labels=c("OL","SG","PW"))+
  scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
  scale_shape_manual(values=c(22,24,22,24,22,24),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))
gg.rich

##water by time
df.div.water <- subset(df.div.exp,type=="Meso. water")
gg.mw.rich.time <- ggplot(df.div.water,aes(x=infusion,y=Observed,fill=inf.temp))+
  geom_boxplot()+
  #geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
  #geom_errorbar(data=df.div.exp.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
  #geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
  #geom_boxplot(outlier.shape=NA,alpha=0.5)+
  facet_wrap(~day)+
  theme_cowplot()+
  ylab("ASV richness")+
  xlab("Infusion")+
  scale_x_discrete(labels=c("OL","SG","PW"))+
  scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
  ggtitle("Mesocosm water")

3.2.2 Simpson’s

df.div.exp.se <- summarySE(df.div.exp,measurevar="InvSimpson",groupvars=c("inf.temp","infusion","temperature","type"))

##dots & error bars
gg.simp <- ggplot(df.div.exp,aes(x=infusion,y=InvSimpson,fill=inf.temp,shape=inf.temp))+
  geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
  geom_errorbar(data=df.div.exp.se,aes(ymax=InvSimpson+se,ymin=InvSimpson-se),width=0.2,color="black",position=position_dodge(width=0.6))+
  geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
  #geom_boxplot(outlier.shape=NA,alpha=0.5)+
  facet_wrap(~type)+
  theme_cowplot()+
  ylab("Simpson's (inv.)")+
  xlab("Infusion")+
  scale_x_discrete(labels=c("OL","SG","PW"))+
  scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
  scale_shape_manual(values=c(22,24,22,24,22,24),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))
gg.simp

gg.mw.simp.time <- ggplot(df.div.water,aes(x=infusion,y=InvSimpson,fill=inf.temp))+
  geom_boxplot()+
  #geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
  #geom_errorbar(data=df.div.exp.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
  #geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
  #geom_boxplot(outlier.shape=NA,alpha=0.5)+
  facet_wrap(~day)+
  theme_cowplot()+
  ylab("Simpson's index (inv.)")+
  xlab("Infusion")+
  scale_x_discrete(labels=c("OL","SG","PW"))+
  scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
  ggtitle("Mesocosm water")

3.2.3 Multi-panel (Figure 2)

ggarrange(gg.rich,gg.simp,common.legend=T,legend="right")

##just water things
ggarrange(gg.mw.rich.time,gg.mw.simp.time,common.legend=T,legend="right",nrow=2)

#ggsave("div.water.time.png",width=6,height=6)

3.3 Statz

3.3.1 Richness

3.3.1.1 Water richness stats

df.div.mw <- subset(df.div.exp,type=="Meso. water")

df.div.mw$day <- as.factor(df.div.mw$day)
str(df.div.mw)
## 'data.frame':    211 obs. of  32 variables:
##  $ sample_name      : chr  "M_D13_OL1_1_S216_L001" "M_D13_OL1_2_S217_L001" "M_D13_OL1_3_S218_L001" "M_D13_OL1_4_S219_L001" ...
##  $ Observed         : num  25 28 22 28 31 31 28 29 27 21 ...
##  $ Shannon          : num  1.63 1.53 1.58 1.53 1.75 ...
##  $ InvSimpson       : num  2.67 2.4 2.44 2.32 2.93 ...
##  $ orgname          : chr  "M_D13_OL1.1" "M_D13_OL1.2" "M_D13_OL1.3" "M_D13_OL1.4" ...
##  $ X16scq           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ totsum           : int  25797 34139 15861 49298 44559 32248 32456 39229 31583 23414 ...
##  $ rspsum           : int  23466 31716 14492 45548 39491 29021 29261 35577 29041 21719 ...
##  $ mesocosm         : chr  "OL1.1" "OL1.2" "OL1.3" "OL1.4" ...
##  $ type             : Factor w/ 2 levels "Mosquitoes","Meso. water": 2 2 2 2 2 2 2 2 2 2 ...
##  $ stage            : chr  "Water" "Water" "Water" "Water" ...
##  $ sex              : chr  "" "" "" "" ...
##  $ date.collected   : chr  "8/27/19" "8/27/19" "8/27/19" "8/27/19" ...
##  $ day              : Factor w/ 3 levels "4","12","20": 2 2 2 2 2 2 2 2 2 2 ...
##  $ abdomen.length.mm: num  NA NA NA NA NA NA NA NA NA NA ...
##  $ wing.length      : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ hindleg.length   : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ special          : chr  "" "" "" "" ...
##  $ dispersal        : chr  "N" "N" "D" "D" ...
##  $ temperature      : chr  "C" "C" "C" "C" ...
##  $ infusion         : chr  "OL" "OL" "OL" "OL" ...
##  $ Wolb_ct1         : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Wolb_ct2         : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Act_ct           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ notes            : chr  "" "" "" "" ...
##  $ newday           : chr  "early" "early" "early" "early" ...
##  $ lib_size         : num  25797 34139 15861 49298 44559 ...
##  $ is.neg           : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
##  $ lib_size_clean   : num  25797 34139 15861 49298 44559 ...
##  $ even             : num  0.506 0.458 0.511 0.46 0.51 ...
##  $ lib_size_trim    : num  25797 34134 15861 49290 44524 ...
##  $ inf.temp         : chr  "OLC" "OLC" "OLC" "OLC" ...
#hist(df.div.mw$Observed)
#shapiro.test(log(df.div.mw$Observed))
shapiro.test(df.div.mw$Observed)
## 
##  Shapiro-Wilk normality test
## 
## data:  df.div.mw$Observed
## W = 0.9775, p-value = 0.001855
#wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),family="poisson",data=df.div.mw)
#wrich2<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),family="compois",data=df.div.mw)
wrich3<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)

#AICtab(wrich1,wrich2,wrich3)

##full model
wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
##no dispersal
wrich1nd<-glmmTMB(Observed~day+temperature+infusion+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
##minus infusion
wrich1ni<-glmmTMB(Observed~day+temperature+dispersal+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
##minus temperature 
wrich1nt<-glmmTMB(Observed~day+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
##minus time 
wrich1nnd<-glmmTMB(Observed~temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)

Anova(wrich1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##                Chisq Df Pr(>Chisq)    
## day          12.9404  2   0.001549 ** 
## temperature   0.6742  1   0.411577    
## infusion    151.4172  2  < 2.2e-16 ***
## dispersal     0.7306  1   0.392704    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Response: Observed
#                Chisq Df Pr(>Chisq)    
# day          13.2024  2   0.001359 ** 
# temperature   0.6742  1   0.411599    
# infusion    152.9448  2  < 2.2e-16 ***
# dispersal     0.7326  1   0.392033    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

AICtab(wrich1,wrich1nd,wrich1ni,wrich1nt,wrich1nnd) #no temp best... almost tied with no dispersal, then full
##           dAIC df
## wrich1nt   0.0 8 
## wrich1nd   0.1 8 
## wrich1     1.3 9 
## wrich1nnd  9.9 7 
## wrich1ni  82.7 7
##interaction?
wrich1.allint <- glmmTMB(Observed~day*temperature*infusion+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.allint1 <- glmmTMB(Observed~day*temperature+infusion+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.allint2 <- glmmTMB(Observed~day+temperature*infusion+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.allint3 <- glmmTMB(Observed~day*infusion+temperature+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw) #second lowest AIC
wrich1.intred <- glmmTMB(Observed~day*infusion+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw) #lowest AIC
wrich1.noint1 <- glmmTMB(Observed~day+infusion+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
wrich1.noint2 <- glmmTMB(Observed~day+infusion+temperature+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)

AICtab(wrich1.allint,wrich1.allint1,wrich1.allint2,wrich1.allint3,wrich1.intred,wrich1.noint1,wrich1.noint2)
##                dAIC df
## wrich1.intred   0.0 11
## wrich1.allint3  2.5 13
## wrich1.allint   5.3 20
## wrich1.noint1  41.8 7 
## wrich1.noint2  43.1 8 
## wrich1.allint2 43.7 10
## wrich1.allint1 45.5 10
#anova(wrich1.int,wrich1.full) #sig 0.001
Anova(wrich1.allint)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##                             Chisq Df Pr(>Chisq)    
## day                       18.4193  2  0.0001001 ***
## temperature                0.8225  1  0.3644497    
## infusion                 177.4337  2  < 2.2e-16 ***
## day:temperature            2.4682  2  0.2911031    
## day:infusion              62.6715  4  7.957e-13 ***
## temperature:infusion       4.1045  2  0.1284452    
## day:temperature:infusion   5.6801  4  0.2243444    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Analysis of Deviance Table (Type II Wald chisquare tests)
# 
# Response: Observed
#                             Chisq Df Pr(>Chisq)    
# day                       18.6821  2  8.775e-05 ***
# temperature                0.8238  1     0.3641    
# infusion                 180.0249  2  < 2.2e-16 ***
# day:temperature            2.3143  2     0.3144    
# day:infusion              62.0786  4  1.060e-12 ***
# temperature:infusion       4.1335  2     0.1266    
# day:temperature:infusion   5.5591  4     0.2346    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Anova(wrich1.intred) #all 0.001, no temperature included
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##                Chisq Df Pr(>Chisq)    
## day           17.390  2  0.0001674 ***
## infusion     164.303  2  < 2.2e-16 ***
## day:infusion  59.023  4  4.653e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(wrich1.allint3) #all 0.001, except temp is 0.3796
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##                 Chisq Df Pr(>Chisq)    
## day           17.4811  2    0.00016 ***
## infusion     167.8537  2  < 2.2e-16 ***
## temperature    0.7595  1    0.38348    
## dispersal      0.7627  1    0.38248    
## day:infusion  58.9948  4  4.717e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#tapply(df.div.mw$log(Observed), df.div.mw$infusion, mean)
#tapply(df.div.mw$log(Observed), df.div.mw$day, mean)

##model checking
shapiro.test(residuals(wrich1.intred)) #awesome
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(wrich1.intred)
## W = 0.9921, p-value = 0.3128
qqnorm(residuals(wrich1.intred))
qqline(residuals(wrich1.intred), col="red")

#plotResiduals(wrich1.int)
wrich1.intred.resid <- simulateResiduals(fittedModel = wrich1.intred, plot = T)

#plot(wrich1.int.resid)
plotResiduals(wrich1.intred.resid, form = df.div.mw$infusion)
## Warning in ensurePredictor(simulationOutput, form): DHARMa:::ensurePredictor:
## character string was provided as predictor. DHARMa has converted to factor
## automatically. To remove this warning, please convert to factor before
## attempting to plot with DHARMa.

plotResiduals(wrich1.intred.resid, form = df.div.mw$day)

##posthoc things
##just without dispersal because not plotting it
##also not plotting day but it should be in the model somewhere already
wrich1.allint3.em <- emmeans(wrich1.allint3,~infusion*temperature)

multcomp::cld(wrich1.allint3.em)
##  infusion temperature emmean    SE  df lower.CL upper.CL .group
##  SW       C             19.1 0.351 198     18.5     19.8  1    
##  SW       H             19.5 0.353 198     18.8     20.1  1    
##  SG       C             21.2 0.355 198     20.5     21.9   2   
##  SG       H             21.5 0.357 198     20.8     22.2   2   
##  OL       C             24.7 0.351 198     24.0     25.4    3  
##  OL       H             25.0 0.353 198     24.3     25.7    3  
## 
## Results are averaged over the levels of: day, dispersal 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
 # infusion temperature emmean    SE  df lower.CL upper.CL .group
 # SW       C             19.1 0.352 199     18.4     19.8  1    
 # SW       H             19.4 0.355 199     18.7     20.1  1    
 # SG       C             21.2 0.356 199     20.5     21.9   2   
 # SG       H             21.5 0.358 199     20.8     22.2   2   
 # OL       C             24.7 0.352 199     24.0     25.4    3  
 # OL       H             25.0 0.355 199     24.3     25.7    3  

3.3.1.2 Water richness stats - rarefied

Same as above but without offset for library size

df.div.mw <- subset(df.div.exp,type=="Meso. water")

df.div.mw$day <- as.factor(df.div.mw$day)
str(df.div.mw)
## 'data.frame':    211 obs. of  32 variables:
##  $ sample_name      : chr  "M_D13_OL1_1_S216_L001" "M_D13_OL1_2_S217_L001" "M_D13_OL1_3_S218_L001" "M_D13_OL1_4_S219_L001" ...
##  $ Observed         : num  25 28 22 28 31 31 28 29 27 21 ...
##  $ Shannon          : num  1.63 1.53 1.58 1.53 1.75 ...
##  $ InvSimpson       : num  2.67 2.4 2.44 2.32 2.93 ...
##  $ orgname          : chr  "M_D13_OL1.1" "M_D13_OL1.2" "M_D13_OL1.3" "M_D13_OL1.4" ...
##  $ X16scq           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ totsum           : int  25797 34139 15861 49298 44559 32248 32456 39229 31583 23414 ...
##  $ rspsum           : int  23466 31716 14492 45548 39491 29021 29261 35577 29041 21719 ...
##  $ mesocosm         : chr  "OL1.1" "OL1.2" "OL1.3" "OL1.4" ...
##  $ type             : Factor w/ 2 levels "Mosquitoes","Meso. water": 2 2 2 2 2 2 2 2 2 2 ...
##  $ stage            : chr  "Water" "Water" "Water" "Water" ...
##  $ sex              : chr  "" "" "" "" ...
##  $ date.collected   : chr  "8/27/19" "8/27/19" "8/27/19" "8/27/19" ...
##  $ day              : Factor w/ 3 levels "4","12","20": 2 2 2 2 2 2 2 2 2 2 ...
##  $ abdomen.length.mm: num  NA NA NA NA NA NA NA NA NA NA ...
##  $ wing.length      : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ hindleg.length   : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ special          : chr  "" "" "" "" ...
##  $ dispersal        : chr  "N" "N" "D" "D" ...
##  $ temperature      : chr  "C" "C" "C" "C" ...
##  $ infusion         : chr  "OL" "OL" "OL" "OL" ...
##  $ Wolb_ct1         : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Wolb_ct2         : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Act_ct           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ notes            : chr  "" "" "" "" ...
##  $ newday           : chr  "early" "early" "early" "early" ...
##  $ lib_size         : num  25797 34139 15861 49298 44559 ...
##  $ is.neg           : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
##  $ lib_size_clean   : num  25797 34139 15861 49298 44559 ...
##  $ even             : num  0.506 0.458 0.511 0.46 0.51 ...
##  $ lib_size_trim    : num  25797 34134 15861 49290 44524 ...
##  $ inf.temp         : chr  "OLC" "OLC" "OLC" "OLC" ...
#hist(df.div.mw$Observed)
#shapiro.test(log(df.div.mw$Observed))
shapiro.test(df.div.mw$Observed)
## 
##  Shapiro-Wilk normality test
## 
## data:  df.div.mw$Observed
## W = 0.9775, p-value = 0.001855
#wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),family="poisson",data=df.div.mw)
#wrich2<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),family="compois",data=df.div.mw)
wrich3<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)

#AICtab(wrich1,wrich2,wrich3)

##full model
wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)
##no dispersal
wrich1nd<-glmmTMB(Observed~day+temperature+infusion+(1|mesocosm), data=df.div.mw)
##minus infusion
wrich1ni<-glmmTMB(Observed~day+temperature+dispersal+(1|mesocosm), data=df.div.mw)
##minus temperature 
wrich1nt<-glmmTMB(Observed~day+infusion+dispersal+(1|mesocosm), data=df.div.mw)
##minus time 
wrich1nnd<-glmmTMB(Observed~temperature+infusion+dispersal+(1|mesocosm), data=df.div.mw)

Anova(wrich1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##                Chisq Df Pr(>Chisq)    
## day           9.4355  2   0.008935 ** 
## temperature   0.4517  1   0.501514    
## infusion    124.7413  2  < 2.2e-16 ***
## dispersal     0.6137  1   0.433416    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Analysis of Deviance Table (Type II Wald chisquare tests)
# 
# Response: Observed
#                Chisq Df Pr(>Chisq)    
# day          10.2547  2   0.005932 ** 
# temperature   0.5020  1   0.478604    
# infusion    133.7700  2  < 2.2e-16 ***
# dispersal     0.5322  1   0.465695    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##same p-values as not rarefied

AICtab(wrich1,wrich1nd,wrich1ni,wrich1nt,wrich1nnd) #no temp best... tied with no dispersal, then full
##           dAIC df
## wrich1nt   0.0 8 
## wrich1nd   0.2 8 
## wrich1     1.5 9 
## wrich1nnd  6.8 7 
## wrich1ni  74.9 7
##interaction?
wrich1.allint <- glmmTMB(Observed~day*temperature*infusion+(1|mesocosm),data=df.div.mw)
wrich1.allint1 <- glmmTMB(Observed~day*temperature+infusion+(1|mesocosm),data=df.div.mw)
wrich1.allint2 <- glmmTMB(Observed~day+temperature*infusion+(1|mesocosm),data=df.div.mw)
wrich1.allint3 <- glmmTMB(Observed~day*infusion+temperature+dispersal+(1|mesocosm),data=df.div.mw) #second lowest AIC
wrich1.intred <- glmmTMB(Observed~day*infusion+(1|mesocosm), data=df.div.mw) #lowest AIC
wrich1.noint1 <- glmmTMB(Observed~day+infusion+(1|mesocosm), data=df.div.mw)
wrich1.noint2 <- glmmTMB(Observed~day+infusion+temperature+(1|mesocosm), data=df.div.mw)

AICtab(wrich1.allint,wrich1.allint1,wrich1.allint2,wrich1.allint3,wrich1.intred,wrich1.noint1,wrich1.noint2) #intred best, same as not rarfied
##                dAIC df
## wrich1.intred   0.0 11
## wrich1.allint3  2.8 13
## wrich1.allint   5.0 20
## wrich1.noint1  45.5 7 
## wrich1.noint2  47.0 8 
## wrich1.allint2 47.4 10
## wrich1.allint1 49.8 10
#anova(wrich1.int,wrich1.full) #sig 0.001
Anova(wrich1.allint)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##                             Chisq Df Pr(>Chisq)    
## day                       13.4995  2   0.001171 ** 
## temperature                0.5746  1   0.448450    
## infusion                 153.0372  2  < 2.2e-16 ***
## day:temperature            1.9400  2   0.379090    
## day:infusion              67.7062  4  6.919e-14 ***
## temperature:infusion       4.5128  2   0.104729    
## day:temperature:infusion   6.4005  4   0.171168    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Analysis of Deviance Table (Type II Wald chisquare tests)
# 
# Response: Observed
#                             Chisq Df Pr(>Chisq)    
# day                       14.7390  2  0.0006302 ***
# temperature                0.6319  1  0.4266476    
# infusion                 164.6677  2  < 2.2e-16 ***
# day:temperature            2.1143  2  0.3474442    
# day:infusion              68.9507  4   3.78e-14 ***
# temperature:infusion       4.6235  2  0.0990855 .  
# day:temperature:infusion   6.1549  4  0.1878773    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Anova(wrich1.intred) #all 0.001, no temperature included
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##                Chisq Df Pr(>Chisq)    
## day           12.743  2   0.001709 ** 
## infusion     141.367  2  < 2.2e-16 ***
## day:infusion  63.747  4  4.725e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(wrich1.allint3) #all 0.001, except temp is 0.44
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##                 Chisq Df Pr(>Chisq)    
## day           12.8152  2   0.001649 ** 
## infusion     143.7560  2  < 2.2e-16 ***
## temperature    0.5250  1   0.468715    
## dispersal      0.6716  1   0.412511    
## day:infusion  63.7223  4  4.781e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#tapply(df.div.mw$log(Observed), df.div.mw$infusion, mean)
#tapply(df.div.mw$log(Observed), df.div.mw$day, mean)

##model checking
shapiro.test(residuals(wrich1.intred)) #awesome
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(wrich1.intred)
## W = 0.99008, p-value = 0.1556
qqnorm(residuals(wrich1.intred))
qqline(residuals(wrich1.intred), col="red")

#plotResiduals(wrich1.int)
wrich1.intred.resid <- simulateResiduals(fittedModel = wrich1.intred, plot = T)

#plot(wrich1.int.resid)
plotResiduals(wrich1.intred.resid, form = df.div.mw$infusion)
## Warning in ensurePredictor(simulationOutput, form): DHARMa:::ensurePredictor:
## character string was provided as predictor. DHARMa has converted to factor
## automatically. To remove this warning, please convert to factor before
## attempting to plot with DHARMa.

plotResiduals(wrich1.intred.resid, form = df.div.mw$day)

##posthoc things
##just without dispersal because not plotting it
##also not plotting day but it should be in the model somewhere already
wrich1.allint3.em <- emmeans(wrich1.allint3,~infusion*temperature)

multcomp::cld(wrich1.allint3.em)
##  infusion temperature emmean    SE  df lower.CL upper.CL .group
##  SW       C             19.2 0.380 198     18.4     19.9  1    
##  SW       H             19.5 0.382 198     18.7     20.2  12   
##  SG       C             21.0 0.383 198     20.2     21.7   23  
##  SG       H             21.2 0.386 198     20.5     22.0    3  
##  OL       C             24.7 0.380 198     23.9     25.4     4 
##  OL       H             25.0 0.382 198     24.2     25.7     4 
## 
## Results are averaged over the levels of: day, dispersal 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
 # infusion temperature emmean    SE  df lower.CL upper.CL .group
 # SW       C             19.1 0.371 198     18.3     19.8  1    
 # SW       H             19.4 0.373 198     18.6     20.1  12   
 # SG       C             20.9 0.374 198     20.1     21.6   23  
 # SG       H             21.1 0.377 198     20.4     21.9    3  
 # OL       C             24.6 0.371 198     23.9     25.4     4 
 # OL       H             24.9 0.373 198     24.2     25.7     4 

3.3.1.3 Mosquito richness stats

df.div.mq <- subset(df.div.exp,type=="Mosquitoes")

hist(df.div.mq$Observed)

shapiro.test(df.div.mq$Observed)
## 
##  Shapiro-Wilk normality test
## 
## data:  df.div.mq$Observed
## W = 0.9207, p-value = 9.183e-09
##different stat families
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
mrich2 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="compois",data=df.div.mq)
mrich3 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="genpois",data=df.div.mq)
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
mrich4 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="poisson",data=df.div.mq)
#mrich4 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="nbinom1",data=df.div.mq)
#mrich5 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="nbinom2",data=df.div.mq)

AICtab(mrich1,mrich2,mrich3,mrich4) #gaussian again
##        dAIC df
## mrich1  0.0 10
## mrich2 62.3 10
## mrich3 63.9 10
## mrich4 64.6 9
##full model
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##no dispersal 
mrich1nd <- glmmTMB(Observed~newday+temperature+infusion+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##minus infusion 
mrich1ni <- glmmTMB(Observed~newday+temperature+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##minus temperature 
mrich1nt <- glmmTMB(Observed~newday+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##minus time 
mrich1nnd <- glmmTMB(Observed~temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mq)
##minus sex
mrich1ns <- glmmTMB(Observed~newday+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)

Anova(mrich1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##              Chisq Df Pr(>Chisq)
## newday      2.8486  2     0.2407
## temperature 0.0704  1     0.7908
## infusion    0.6272  2     0.7308
## dispersal   0.9250  1     0.3362
## sex         0.0527  1     0.8185
shapiro.test(residuals(mrich1))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mrich1)
## W = 0.93556, p-value = 1.302e-07
qqnorm(residuals(mrich1))
qqline(residuals(mrich1), col="red")

mrich1.resid <- simulateResiduals(fittedModel = mrich1, plot = T) 

##posthoc things
mrich1.em <- emmeans(mrich1,~infusion*temperature)
multcomp::cld(mrich1.em)
##  infusion temperature emmean    SE  df lower.CL upper.CL .group
##  OL       H             8.05 0.395 185     7.27     8.83  1    
##  OL       C             8.17 0.323 185     7.53     8.80  1    
##  SG       H             8.19 0.425 185     7.35     9.03  1    
##  SG       C             8.31 0.465 185     7.39     9.23  1    
##  SW       H             8.51 0.454 185     7.61     9.40  1    
##  SW       C             8.63 0.531 185     7.58     9.67  1    
## 
## Results are averaged over the levels of: newday, dispersal, sex 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
##no differences

3.3.1.4 Mosquito richness stats - rarefied

Same as above but without offset for library size

df.div.mq <- subset(df.div.exp,type=="Mosquitoes")

hist(df.div.mq$Observed)

shapiro.test(df.div.mq$Observed)
## 
##  Shapiro-Wilk normality test
## 
## data:  df.div.mq$Observed
## W = 0.9207, p-value = 9.183e-09
##different stat families
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),data=df.div.mq)
mrich2 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="compois",data=df.div.mq)
mrich3 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="genpois",data=df.div.mq)
mrich4 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="poisson",data=df.div.mq)
#mrich4 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="nbinom1",data=df.div.mq)
#mrich5 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="nbinom2",data=df.div.mq)

AICtab(mrich1,mrich2,mrich3,mrich4) #compois 
##        dAIC df
## mrich2  0.0 10
## mrich3  6.5 10
## mrich4 18.3 9 
## mrich1 18.9 10
##full model
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="compois",data=df.div.mq)
##no dispersal 
mrich1nd <- glmmTMB(Observed~newday+temperature+infusion+sex+(1|mesocosm),family="compois",data=df.div.mq)
##minus infusion 
mrich1ni <- glmmTMB(Observed~newday+temperature+dispersal+sex+(1|mesocosm),family="compois",data=df.div.mq)
##minus temperature 
mrich1nt <- glmmTMB(Observed~newday+infusion+dispersal+sex+(1|mesocosm),family="compois",data=df.div.mq)
##minus time 
mrich1nnd <- glmmTMB(Observed~temperature+infusion+dispersal+sex+(1|mesocosm),family="compois", data=df.div.mq)
##minus sex
mrich1ns <- glmmTMB(Observed~newday+temperature+infusion+dispersal+(1|mesocosm),family="compois",data=df.div.mq)

Anova(mrich1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##              Chisq Df Pr(>Chisq)
## newday      3.9718  2     0.1373
## temperature 0.0016  1     0.9677
## infusion    0.2076  2     0.9014
## dispersal   0.6153  1     0.4328
## sex         0.0050  1     0.9435
shapiro.test(residuals(mrich1))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mrich1)
## W = 0.93647, p-value = 1.55e-07
qqnorm(residuals(mrich1))
qqline(residuals(mrich1), col="red")

mrich1.resid <- simulateResiduals(fittedModel = mrich1, plot = T) 

##posthoc things
mrich1.em <- emmeans(mrich1,~infusion*temperature)
multcomp::cld(mrich1.em)
##  infusion temperature emmean     SE  df asymp.LCL asymp.UCL .group
##  OL       C             2.09 0.0428 Inf      2.00      2.17  1    
##  OL       H             2.09 0.0504 Inf      1.99      2.19  1    
##  SG       C             2.11 0.0606 Inf      1.99      2.23  1    
##  SG       H             2.11 0.0558 Inf      2.00      2.22  1    
##  SW       C             2.11 0.0653 Inf      1.99      2.24  1    
##  SW       H             2.12 0.0573 Inf      2.00      2.23  1    
## 
## Results are averaged over the levels of: newday, dispersal, sex 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
##no differences

3.3.2 Simpson’s

3.3.2.1 Water Simpson’s stats

hist(df.div.mw$InvSimpson)

shapiro.test(df.div.mw$InvSimpson)
## 
##  Shapiro-Wilk normality test
## 
## data:  df.div.mw$InvSimpson
## W = 0.96756, p-value = 9.003e-05
hist(log(df.div.mw$InvSimpson))

shapiro.test(log(df.div.mw$InvSimpson))
## 
##  Shapiro-Wilk normality test
## 
## data:  log(df.div.mw$InvSimpson)
## W = 0.95673, p-value = 5.166e-06
##normalizing because residuals bad below
bestNormalize(df.div.mw$InvSimpson) #orderNorm
## Best Normalizing transformation with 211 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 1.5829
##  - Box-Cox: 1.3786
##  - Center+scale: 1.3483
##  - Double Reversed Log_b(x+a): 1.498
##  - Exp(x): 3.4691
##  - Log_b(x+a): 1.6639
##  - orderNorm (ORQ): 1.1158
##  - sqrt(x + a): 1.3829
##  - Yeo-Johnson: 1.385
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## orderNorm Transformation with 211 nonmissing obs and no ties 
##  - Original quantiles:
##    0%   25%   50%   75%  100% 
## 1.359 1.947 2.669 3.225 5.130
df.div.mw$simp.norm <- bestNormalize(df.div.mw$InvSimpson)$x.t
shapiro.test(df.div.mw$simp.norm)
## 
##  Shapiro-Wilk normality test
## 
## data:  df.div.mw$simp.norm
## W = 0.96756, p-value = 9.003e-05
##just replaced InvSimpson with simp.norm
##full model
wsimp1 <- glmmTMB(simp.norm~day+temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)
##no dispersal
wsimp1nd <- glmmTMB(simp.norm~day+temperature+infusion+(1|mesocosm),data=df.div.mw)
##minus infusion
wsimp1ni <- glmmTMB(simp.norm~day+temperature+dispersal+(1|mesocosm),data=df.div.mw)
##minus temperature 
wsimp1nt <- glmmTMB(simp.norm~day+infusion+dispersal+(1|mesocosm),data=df.div.mw)
##minus time 
wsimp1nnd <- glmmTMB(simp.norm~temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)

Anova(wsimp1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: simp.norm
##                Chisq Df Pr(>Chisq)    
## day          45.4454  2  1.354e-10 ***
## temperature  48.9643  1  2.607e-12 ***
## infusion    276.6117  2  < 2.2e-16 ***
## dispersal     0.0291  1     0.8644    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Analysis of Deviance Table (Type II Wald chisquare tests)
# 
# Response: simp.norm
#                Chisq Df Pr(>Chisq)    
# day          49.7901  2  1.542e-11 ***
# temperature  41.3001  1  1.306e-10 ***
# infusion    229.1055  2  < 2.2e-16 ***
# dispersal     0.0091  1      0.924    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##~same results with normalized results
##same results with rarefied

AICtab(wsimp1,wsimp1nd,wsimp1ni,wsimp1nt,wsimp1nnd)
##           dAIC  df
## wsimp1nd    0.0 8 
## wsimp1      2.0 9 
## wsimp1nt   37.3 8 
## wsimp1nnd  37.5 7 
## wsimp1ni  112.4 7
##no dispersal good, then full model

##interaction?
wsimp1.int1 <- glmmTMB(simp.norm~day*infusion*temperature+(1|mesocosm), data=df.div.mw)
wsimp1.int2 <- glmmTMB(simp.norm~day*infusion+temperature+(1|mesocosm), data=df.div.mw) #best AIC
wsimp1.int3 <- glmmTMB(simp.norm~day+infusion*temperature+(1|mesocosm), data=df.div.mw)

AICtab(wsimp1.int1,wsimp1.int2,wsimp1.int3,wsimp1nd) #int better
##             dAIC df
## wsimp1.int2  0.0 12
## wsimp1.int1  2.1 20
## wsimp1.int3 33.6 10
## wsimp1nd    36.8 8
#interaction with day & infusion

Anova(wsimp1.int2) #all 0.001 sig level, same with rare
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: simp.norm
##                Chisq Df Pr(>Chisq)    
## day           62.772  2  2.340e-14 ***
## infusion     274.807  2  < 2.2e-16 ***
## temperature   47.216  1  6.358e-12 ***
## day:infusion  53.210  4  7.703e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Response: simp.norm
#                Chisq Df Pr(>Chisq)    
# day           70.733  2  4.370e-16 ***
# infusion     227.491  2  < 2.2e-16 ***
# temperature   39.606  1  3.107e-10 ***
# day:infusion  58.541  4  5.874e-12 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

##all of the below was bad before normalizing
shapiro.test(residuals(wsimp1))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(wsimp1)
## W = 0.95063, p-value = 1.211e-06
qqnorm(residuals(wsimp1))
qqline(residuals(wsimp1), col="red")

wsimp1.resid <- simulateResiduals(fittedModel = wsimp1, plot = T) 

##posthoc things
wsimp1.em <- emmeans(wsimp1,~infusion*temperature)
multcomp::cld(wsimp1.em)
##  infusion temperature emmean     SE  df lower.CL upper.CL .group
##  SG       C           -1.298 0.0882 202   -1.472  -1.1239  1    
##  SG       H           -0.682 0.0887 202   -0.857  -0.5074   2   
##  OL       C           -0.109 0.0875 202   -0.282   0.0635    3  
##  SW       C            0.464 0.0875 202    0.292   0.6368     4 
##  OL       H            0.507 0.0880 202    0.333   0.6801     4 
##  SW       H            1.080 0.0880 202    0.906   1.2533      5
## 
## Results are averaged over the levels of: day, dispersal 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
 # infusion temperature emmean     SE  df lower.CL upper.CL .group
 # SG       C           -1.255 0.0934 202   -1.439  -1.0708  1    
 # SG       H           -0.656 0.0939 202   -0.841  -0.4712   2   
 # OL       C           -0.104 0.0927 202   -0.287   0.0786    3  
 # SW       C            0.442 0.0927 202    0.259   0.6244     4 
 # OL       H            0.495 0.0932 202    0.311   0.6783     4 
 # SW       H            1.040 0.0932 202    0.857   1.2240      5

3.3.2.2 Mosquito Simpson’s

From manuscript pre-revisions: “In contrast, the Simpson’s index, which integrates evenness into a measure of alpha diversity, of the adult mosquito microbiome was lower in male mosquitoes relative to female mosquitoes (p<0.0001). This trend was associated with the increased dominance of wAlbB in males. The Simpson’s index of the adult mosquito microbiome did not vary significantly with the dispersal, aquatic chemistry, temperature, or time of emergence.”

hist(df.div.mq$InvSimpson)

df.div.mq$simp.norm <- bestNormalize(df.div.mq$InvSimpson)$x.t #orderNorm
hist(df.div.mq$simp.norm)

##full model
msimp.all <- glmmTMB(simp.norm~newday+temperature+infusion+dispersal+sex+(1|mesocosm), data=df.div.mq)
##without sex
msimp.nosex <- glmmTMB(simp.norm~day+temperature+infusion+dispersal+(1|mesocosm), data=df.div.mq)
##without dispersal
msimp.nodis <- glmmTMB(simp.norm~newday+temperature+infusion+sex+(1|mesocosm), data=df.div.mq)
##without infusion
msimp.noinf <- glmmTMB(simp.norm~newday+temperature+dispersal+sex+(1|mesocosm), data=df.div.mq)
##without temperature
msimp.notem <- glmmTMB(simp.norm~newday+infusion+dispersal+sex+(1|mesocosm), data=df.div.mq)
##without time
msimp.notim <- glmmTMB(simp.norm~temperature+infusion+dispersal+sex+(1|mesocosm), data=df.div.mq)

Anova(msimp.all)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: simp.norm
##               Chisq Df Pr(>Chisq)    
## newday       2.3243  2    0.31281    
## temperature  0.2966  1    0.58605    
## infusion     2.7471  2    0.25321    
## dispersal    2.8196  1    0.09312 .  
## sex         67.6024  1    < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Analysis of Deviance Table (Type II Wald chisquare tests)
# 
# Response: simp.norm
#               Chisq Df Pr(>Chisq)    
# newday       2.2911  2    0.31805    
# temperature  0.2987  1    0.58469    
# infusion     2.7011  2    0.25909    
# dispersal    2.8332  1    0.09234 .  
# sex         67.5319  1    < 2e-16 ***

AICtab(msimp.all,msimp.nosex,msimp.nodis,msimp.noinf,msimp.notem,msimp.notim)
##             dAIC df
## msimp.notem  0.0 9 
## msimp.notim  0.0 8 
## msimp.noinf  0.4 8 
## msimp.all    1.7 10
## msimp.nodis  2.4 9 
## msimp.nosex 59.1 8
#no time & no temp tied for best, then no infusion

##all of the below was bad before normalizing
shapiro.test(residuals(msimp.all))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(msimp.all)
## W = 0.98795, p-value = 0.09731
qqnorm(residuals(msimp.all))
qqline(residuals(msimp.all), col="red")

msimp.all.resid <- simulateResiduals(fittedModel = msimp.all, plot = T) 

##posthoc things
msimp.all.em <- emmeans(msimp.all,~infusion*temperature)
multcomp::cld(msimp.all.em)
##  infusion temperature  emmean    SE  df lower.CL upper.CL .group
##  SG       H           -0.2061 0.152 185   -0.506   0.0939  1    
##  SG       C           -0.1164 0.167 185   -0.446   0.2136  1    
##  OL       H           -0.0764 0.143 185   -0.359   0.2063  1    
##  OL       C            0.0132 0.115 185   -0.215   0.2409  1    
##  SW       H            0.1483 0.162 185   -0.172   0.4682  1    
##  SW       C            0.2379 0.193 185   -0.142   0.6182  1    
## 
## Results are averaged over the levels of: newday, dispersal, sex 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
#no diffs

3.4 Infusion water things

df.div.inf <- subset(df.div,type=="Infusion water")

df.div.inf.se <- summarySE(df.div.inf,measurevar="Observed",groupvars=c("infusion"))
df.div.inf.se
##   infusion N  Observed        sd        se        ci
## 1       OL 3 44.666667 4.1633320 2.4037009 10.342290
## 2       SG 3 20.666667 2.0816660 1.2018504  5.171145
## 3       SW 3  3.333333 0.5773503 0.3333333  1.434218
# gg.iw.rich <- ggplot(df.div.inf,aes(x=infusion,y=Observed,fill=infusion))+
#   geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
#   geom_errorbar(data=df.div.inf.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
#   geom_point(data=df.div.inf.se,size=2.5,position=position_dodge(width=0.6))+
#   #geom_boxplot(outlier.shape=NA,alpha=0.5)+
#   theme_cowplot()+
#   ggtitle("Infusion water")
# gg.iw.rich

4 Diversity correlations

4.1 Diversity x day

Doesn’t look interesting

ggplot(df.div.mq,aes(x=day,y=Observed))+
  geom_point()+
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

ggplot(df.div.mq,aes(x=day,y=InvSimpson))+
  geom_point()+
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

4.2 Diversity - mesocosm vs. mosquito

library(reshape)
## 
## Attaching package: 'reshape'
## The following object is masked from 'package:Matrix':
## 
##     expand
## The following object is masked from 'package:cowplot':
## 
##     stamp
## The following objects are masked from 'package:plyr':
## 
##     rename, round_any
meso.corr <- merge(df.div.mq,df.div.mw,by="mesocosm")
ggplot(meso.corr,aes(x=InvSimpson.x,y=InvSimpson.y))+
  geom_point()+
  #facet_wrap(~day.y)+
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

ggplot(meso.corr,aes(x=Observed.x,y=Observed.y))+
  geom_point()+
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'